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THE  COMMINUTION  BEHAVIOR  OF  LIQUIDS 


1.  INTRODUCTION 

The  aim  of  this  investigation  is  to  assess  the  role  of  nonlinear  viscoelasticity  in  the 
process  of  liquid  comminution  under  explosive  stress.  Ultimately  the  goal  of  the  program  is  to 
evaluate  the  hypothesis  that  the  crucial  factor  in  determining  the  comminution  and  dispersal 
characteristics  of  liquids  under  these  conditions  is  the  nonlinear  dynamical  response  behavior  of 
the  liquids.  In  attacking  the  problem,  it  is  essential  first  to  develop  a  model  of  comminution 
that  incorporates  the  mechanism  of  liquid  response  and  failure  when  subjected  to  explosively 
generated  stresses.  The  central  theme  of  this  report  is  the  proposal  of  such  a  model,  the  ex¬ 
ploration  of  its  consequences,  and,  where  possible,  a  comparison  of  predictions  based  upon  it 
with  experimental  data.  This  model  calculation  forms  the  first  part  of  the  report.  Following 
this  is  a  brief  section  which  isolates  those  elements  of  the  proposed  mechanism  that  are  con¬ 
trolled  by  the  nonlinear  viscoelastic  characteristics  of  the  liquid.  Lastly  we  address  the  issue  of 
utilizing  computer-simulation  molecular-dynamics  (MD)  experiments  as  a  tool  for  probing  these 
nonlinear  characteristics;  we  conclude  that  MD  experiments  offer  significant  potential  for 
contributing  to  the  understanding  and  control  of  the  processes  of  explosive  liquid  comminution 
and  dispersal. 

2.  BACKGROUND 

The  experimental  study  that  provides  the  principal  impetus  and  guidance  for  addressing 
comminution  behavior  is  that  of  Gerber  and  Stuempfle.*  In  this  work  a  high-explosive  projector 
was  employed  to  produce  dissemination  of  a  liquid  cloud  around  a  preferential  vertical  axis.  A 
three-dimensional  sampling  array  was  used  to  measure  the  mass  distribution  of  the  cloud  follow¬ 
ing  its  emergence  from  the  projector.  The  mass  distribution  along  the  vertical  z-axis  was 
described  by  a  Weibull  distribution, 

f(z)  =  (0/z)(z/zq)^  exp[-(z/z0)^J 

such  that  f(z)dz  gives  the  mass  fraction  lost  in  the  height  range  z  to  z  +  dz  above  the 
projector.  In  the  Weibull  form  Zq  is  a  characteristic  length  parameter  and  0  is  a  shape  parame¬ 
ter.  Gerber  and  Stuempfle  present  their  data  in  terms  of  four  parameters  z^q  —  the  median  of 
the  distribution  —  and  the  first  three  moments  of  the  distribution.  In  terms  of  the  parameters 
of  the  distribution 


z50  _z0^7?  ^ 

and  the  mean  or  first  moment  is 

<z>  =  zQr (1/6 )  /g  . 

Values  of  ?  fall  in  the  range  from  about  0.6  to  3  with  most  of  the  more  than  30  liquids  studied 
being  described  by  ? 's  near  2;  <^>  values  varied  from  about  65cm  to  about  335cm,  Generally 
one  finds  that  larger  values  of  <[z)>  are  accompanied  by  larger  values  of  (J.  A  major  theme  in 
what  follows  is  the  attempt  to  understand  those  liquid  properties  that  are  operative  in  the 


*  Gerber,  B.V.  and  Stuempfle,  A.D.  TSD  004.7-1 1 /Vol.  1.  A  New  Experimental  Technique  for 
Studying  the  Explosive  Comminution  of  Liquids.  1976  Army  Science  Conference 
Proceedings.  Vol.  I.  September  1976.  Unclassified  Report. 
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determination  of  the  average  height  <z>.  Exploring  the  details  of  the  mass  distribution  pattern 
is  beyond  the  scope  of  this  work. 


4.  PRELIMINARY  CONSIDERATIONS 

In  this  section  we  give  a  crude  description  of  the  liquid  comminution  process.  The 
principal  result  of  the  section  is  an  equation  giving  the  average  height  <z>  in  terms  of  the  size 
of  the  particles  or  droplets  that  compose  the  cloud. 

We  consider  the  effect  of  the  explosion  on  the  liquid  to  be  roughly  separable  into  two 
major  components: 

(a)  The  explosive  release  of  energy  generates  a  shock  wave  that  in  passing  through 
the  liquid  produces  large  tensile,  compressional,  and  shear  stresses.  These  cause  local  internal 
fissures  within  the  liquid  so  that  it  comminutes  into  many  small  liquid  fragments. 

(b)  The  explosion  produces  a  rapidly  expanding  gas  that  intrudes  itself  into  the 
stress-produced  fissures,  thereby  preventing  their  healing,  and  causing  the  liquid  fragments  to 
be  expelled  from  the  projector  with  a  vertical  velocity  Vq. 

The  dependence  of  the  initial  velocitv  on  liquid  parameters  can  be  estimated  by  the  fol¬ 
lowing  sort  of  rough  argument.  We  shall  assume  that  relatively  little  of  the  energy  of  explosion 
is  used  in  the  formation  of  the  droplets  so  that  most  of  the  energy  E  is  transferred  to  the 
kinetic  energy  of  the  droplets.  The  initial  velocity  imparted  to  the  droplets  then  becomes 


V0  =  (2e/tn)1/2xp'1/2  (1) 

where  m  is  the  mass  of  a  droplet,  G  is  the  kinetic  energy  per  droplet,  and  p  is  the  liquid 
densitv.  Assuming  that  E  is  the  same  for  each  explosion,  it  follows  that  the  initial  velocity  of 
the  droplets  is  dependent  only  on  the  density  of  the  liquid. 

The  dynamical  problem  for  the  projected  droplets  can  now  be  worked  out 
straightforwardly,  at  least  in  lowest  order  approximation.  We  assume  that  the  droplets  are 
spherical  and  that  the  size  distribution  is  characterized  by  a  function  f(a)  which  gives  the  proba¬ 
bility  of  a  droplet's  radius  falling  in  the  range  a  to  a+da.  We  take  the  forces  acting  on  a  droplet 
of  radius  a  to  be 


and 


(a)  gravity:  ^Tra^pg, 

(b)  air  resistance:  -67r*?av. 


In  these  expressions  g  is  the  acceleration  due  to  gravity,  r\  is  the  viscosity  of  the  medium  fair) 
into  which  the  droplet  cloud  is  projected,  and  v  is  the  instantaneous  droplet  velocitv  in  the 
vertical  direction.  Newton's  law  then  takes  the  form 


dv/cU  +  v/x  = -g, 


(2) 


where  for  a  droplet  of  radius  a 


T  =  2pa^/9n. 


R 


(3) 
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The  solution  of  (2)  is  just 


(4) 


v  =  v0  exp(-t/x)-gx  [1  -  exp  (-t/x)]  . 

Integrating  this  again  gives  the  equation  of  motion 

z(t)  =  v0xtl-exp  (-t/x)  +  gx2  1  -  Y-exP  (-t/x)].  (5) 

The  height  to  which  a  droplet  rises,  z,  can  then  be  found  by  setting  v=0  in  (4)  and  solving  for  the 
time: 
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(6) 


Putting  z  =  2  (t  )  gives 

0 

2 

z=v0x-gx  HnO+v0/gx).  (7) 

o 

For  typical  values  of  r? ,  p,  and  a  (v  =L  0.18  mP,Pc^l.3  g/cm  ,  a^25  pm),T  — .01  sec.  Using  3 
x  10*  cm/s ec  and  g  =  980  cm/sec^  one  gets  t  ~  .08  sec,  which  is  of  the  correct  order  of 
magnitude  to  what  is  experimentally  observed.  Moreover  vQ»gr  so  that  approximately 

z  -  vqT.  (8) 

Equation  (8)  is  accurate  to  better  than  a  percent  or  so  for  the  full  range  of  liquids  and  experi¬ 
mental  conditions  studied  by  Gerber  and  Stuempfle.  Using  equation  (3)  and  averaging  gives 

<z>  =  (2vQ/9n)p<a2’>,  (9) 

where  <\  a  y  is  the  mean  square  droplet  radius. 


4.  A  PHYSICAL  MODEL  FOR  COMMINUTION 

The  central  issue  is  clear  from  equation  (9):  What  are  the  processes  -  and  thus  the 
liquid  characteristics  -  that  determine  the  mean  square  droplet  size?  We  assume  that  the 
mechanism  of  droplet  formation  is  mainly  a  mechanical  (as  opposed  to  a  thermal)  one.  At  any 
instant  even  an  equilibrium  liquid  exhibits  local  deviations  from  homogeneity  and  isotropy.  One 
characterization  of  the  deviations  from  homogeneity  is  the  distribution  of  fluctuations  in  the 
density.  From  classical  statistical  mechanics  the  probability  of  occurrence  of  a  density 
fluctuation  in  the  range  Apto  Ap+  dp  is 
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,  2 

where  <  Ap  >  is  the  mean  square  density  fluctuation  given  by 

<  Ap2>  =  p2  (kT/V)/K,  (11) 


where  T  and  p  are  the  temperature  and  average  density  of  the  system,  V  is  the  volume  of  the 
region  under  consideration,  k  is  Boltzmann's  constant,  and  K  is  the  equilibrium  bulk  modulus  of 
the  material. 

We  suggest  that  under  the  influence  of  the  large  explosively  generated  stresses,  these 
local  inhomogeneities  will  behave  as  structural  flaws  in  the  material:  shear  stresses  tend  to 
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elongate  the  inhomogeneities  by  "organizing"  the  liquid  structure*  while  the  tensile  stresses 
cause  the  flaws  to  propagate  at  the  stress  intensification  points.  The  situation  is  depicted 
pictorially  in  figure  1.  However,  in  a  liquid,  even  at  the  elevated  pressures  produced  by  the 
explosive  shock,  one  expects  that  the  rate  at  which  the  flaws  heal  and  the  rate  at  which  they 
propagate  via  growth  at  their  tips  should  be  comparable;  consequently,  one  is  not  dealing  with  a 
classical  brittle  fracture  situation,  but  rather  one  involving  relatively  little  growth.  Thus  the 
maximum  flaw  dimension  should  be  comparable  with  the  size  of  the  particles  that  result  from 
the  explosive  comminution. 

It  is  the  processes  described  in  this  last  paragraph  that  will  depend  sensitively  on  the 
nonlinear  viscoelastic  response  characteristics  of  the  liquid  system.  Because  the  large  stresses 
to  which  the  material  is  subjected  are  "switched  on"  essentially  instantaneously,  the  nature  of 
the  behavior  will  be  governed  in  large  measure  by  the  characteristic  time  scales  of  the  shear 
and  structural  responses  of  the  liquid,  i.e.,  the  shear  and  structural  relaxation  times.  The 
specific  picture  that  we  have  described  above  is  an  oversimplified  one  in  which  we  have  as¬ 
sumed  that  the  .hear  relaxation  time  is  negligibly  small  in  comparison  with  the  structural 
relaxation  time.  Only  in  this  limiting  situation  is  one  able  to  portray  the  flaw  development  and 
liquid  fracture  processes  without  reference  to  the  specific  viscoelastic  and  strength  parameters 
that  describe  the  material.  Generally  speaking  these  extreme  conditions  do  not  obtain;  the 
shear  relaxation  time,  while  smaller  than  the  structural  time,  is  of  the  same  order  of  magni¬ 
tude.  Consequently,  the  sequential  nature  of  the  mechanisms  that  operate  in  comminution  that 
was  assumed  above  —  the  development  of  large  shear  strains,  followed  by  structural  reorganiza¬ 
tion  (flaw  development),  followed  by  liquid  fracture  from  tensile  stresses  —  is  not  nearly  so 
precisely  defineable.  The  various  mechanisms  do  not  proceed  strictly  serially,  but  rather  couple 
in  a  complicated  interactive  and,  perhaps,  competitive  fashion.  Depending  upon  the  relative 
time  scales  and  specific  forms  of  the  various  response  mechanisms,  this  coupling  can  either 
reinforce  the  liquid  comminution  process  (leading  to  smaller  particle  sizes)  or  interfere  with  it 
^resulting  in  larger  particle  sizes).  Nevertheless,  the  elementary  picture  that  we  are  examining 
is  a  useful  one  for  it  provides  a  general  framework  in  which  to  analyze  and  compare  the 
behavior  of  different  liquids  acted  on  by  explosive  stresses,  as  well  as  highlighting  the  crucial 
liquid  properties  that  govern  (at  least  in  lowest  order  approximation)  the  cloud  dispersal  charac¬ 
teristics, 

In  order  to  estimate  the  droplet  size,  we  can  employ  the  following  sort  of  element arv 
notions.  We  consider  that  in  order  to  generate  droplets  of  charactoristir  dimension  =  a,  a  region 
of  the  liquid  of  volume  must  contain  a  flaw  whose  volume  exceeds  some  critical  volume 

V  ,  The  action  of  the  explosive  shear  stresses  is  to  elongate  and  shape  the  flaw  [into  roughly  a 
"pillow"  shape  (see  figure  2)]  so  that  it  can  be  described  by  dimensions  a,  b,  and  c  (a>b>c).  The 
minimum  flaw  volume  V  =  abc,  that  will  lead  to  the  production  of  droplets  with  dimensions  on 
the  order  of  "a"  will  result  if  one  has  a  local  density  fluctuation 

Ap  >=pVc/v  =  A pr  (12) 

in  the  volume  V.  The  liquid  will  comminute  into  particles  of  this  size  only  if  a  fluctuation  of 
such  magnitude  occurs  with  significant  probabilitv.  The  condition  for  this  to  obtain  is 

Apc«<  ApZ>1/2:  (13) 


*  Heyes ~ F)T~M T,  Kim,  J.  J.,  Montrose,  C.  J.  and  Litovitz,  T.  A.  Time-Dependent  Non-Linear 
Shear  Stress  Effects  in  Simple  Liquids:  A  Molecular  Dynamics  Study.  J.  Chem.  Phys.  73, 
3987  (1980). 
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Figure  l.A.  A  Schematic  Representation  of  Liquid  Flaws  Elongated  under  the  Action  of  Large 

Shear  Stresses 


Figure  1. 


•  Liquid  Flaws  That  Have  Grown  and  Intersected  to  Form  Comminution  Droplets 
under  the  Action  of  Tensile  Stresses. 


that  is,  the  critical  density  fluctuation  must  be  comparable  with  the  PMS  density  fluctuation. 
Substituting  equations  (11)  and  (12)  in  equation  (13)  and  using  V«a3  and  V  =  abc  leads  to  an 
approximate  expression  for  the  value  of  a: 


a  =  (be)2  K/kT. 

Recalling  equation  (9),  this  leads  to  an  equation  for  the  mean  height  of  the  droplets: 


(14) 


<z>  =  (2/9n)  vQp  (be)4  (K/kT)2.  (15) 

In  order  that  the  flaws  do  not  heal  themselves  too  rapidly,  it  is  essential  that  the  smallest  flaw 
dimension  c  be  at  least  on  the  order  of  a  molecular  separation  distance;  with  this  in  mind  we 
have 


-1/3 

CaP  (16) 

-1/2 

and  recalling  that  v0  x>  p  7  ,  we  may  collect  together  all  the  density-dependent  quantities  on 
the  right  hand  side  of  equation  (15).  The  resulting  expression  takes  the  form 


<z >  oc  p-5/6  y/2.  (17) 


5.  COMPARISON  WITH  EXPERIMENT 

Equation  (17)  is  in  a  form  that  can  be  compared  with  experimental  results.  The 
comparison  is  shown  in  figure  3  for  six  liquids  at  20°C.  Many  of  the  data  for  this  figure  were 
taken  from  the  work  of  Gerber  and  Stuempfle*  and  are  presented  in  the  table  which  will 
follow.  As  is  apparent  there  is  fairlv  good  agreement  between  these  data  and  the  prediction  of 
equation  (17)  (the  slope  of  the  line  shown  in  figure  3  is  65.5  g^7”  cm'^^GPa"  ').  This  offers 
fairly  strong  support  for  the  essential  correctness  of  the  view  of  the  comminution  mechanism 
that  was  described  above.  It  is  also  clear  that  the  modifications  of  this  overly  simplified  view 
that  are  required  if  one  correctly  treats  the  dynamical  aspects  of  the  liquid  response  are  not 
understood  in  any  kind  of  quantitative  fashion.  Some  of  the  molecular  dynamics  results 
(presented  in  the  appendix)  give  a  fair  indication  of  the  type  of  phenomena  that  must  be 
considered  —  especially  in  regard  to  liquid  structural  changes  occasioned  by  large  shear  stress  — 
but  efforts  to  apply  these  results  quantitatively  to  comminution  are  at  this  point  premature. 

Of  the  six  liquids  which  were  considered  above,  except  for  carbon  tetrachloride, 
comminution  data  have  also  been  obtained  at  130°C.  The  variation  that  one  would  expect  from 
our  model  is  given  by  equation  (15);  the  major  effect  appears  to  be  the  factor  (1/kT)  .  On  the 
basis  of  this  one  would  expect  the  ratio  of  the  mean  height  at  20°C  to  that  at  130°C  would  be 
about  (403/293)^^1 .9.  The  experimentally  determined  ratios  for  the  liquids  tha*  we  have  con¬ 
sidered  fall  into  the  range  of  about  1.0  to  1.5. 


*Gerber,  B.  V.,  and  Stuempfle,  A.  D.,  op.  cit. 


that  is,  the  critical  density  fluctuation  must  be  comparable  with  the  PMS  density  fluctuation. 
Substituting  equations  (11)  and  (12)  in  equation  (13)  and  using  V«a3  and  V  =  abc  leads  to  an 
approximate  expression  for  the  value  of  a: 


a  =  (be)2  K/kT. 

Recalling  equation  (9),  this  leads  to  an  equation  for  the  mean  height  of  the  droplets: 


(14) 


<z>  =  (2/9n)  vQp  (be)4  (K/kT)2.  (15) 

In  order  that  the  flaws  do  not  heal  themselves  too  rapidly,  it  is  essential  that  the  smallest  flaw 
dimension  c  be  at  least  on  the  order  of  a  molecular  separation  distance;  with  this  in  mind  we 
have 


-1/3 

CaP  (16) 

-1/2 

and  recalling  that  v0  x>  p  7  ,  we  may  collect  together  all  the  density-dependent  quantities  on 
the  right  hand  side  of  equation  (15).  The  resulting  expression  takes  the  form 


<z >  oc  p-5/6  y/2.  (17) 


5.  COMPARISON  WITH  EXPERIMENT 

Equation  (17)  is  in  a  form  that  can  be  compared  with  experimental  results.  The 
comparison  is  shown  in  figure  3  for  six  liquids  at  20°C.  Many  of  the  data  for  this  figure  were 
taken  from  the  work  of  Gerber  and  Stuempfle*  and  are  presented  in  the  table  which  will 
follow.  As  is  apparent  there  is  fairlv  good  agreement  between  these  data  and  the  prediction  of 
equation  (17)  (the  slope  of  the  line  shown  in  figure  3  is  65.5  g^7”  cm'^^GPa"  ').  This  offers 
fairly  strong  support  for  the  essential  correctness  of  the  view  of  the  comminution  mechanism 
that  was  described  above.  It  is  also  clear  that  the  modifications  of  this  overly  simplified  view 
that  are  required  if  one  correctly  treats  the  dynamical  aspects  of  the  liquid  response  are  not 
understood  in  any  kind  of  quantitative  fashion.  Some  of  the  molecular  dynamics  results 
(presented  in  the  appendix)  give  a  fair  indication  of  the  type  of  phenomena  that  must  be 
considered  —  especially  in  regard  to  liquid  structural  changes  occasioned  by  large  shear  stress  — 
but  efforts  to  apply  these  results  quantitatively  to  comminution  are  at  this  point  premature. 

Of  the  six  liquids  which  were  considered  above,  except  for  carbon  tetrachloride, 
comminution  data  have  also  been  obtained  at  130°C.  The  variation  that  one  would  expect  from 
our  model  is  given  by  equation  (15);  the  major  effect  appears  to  be  the  factor  (1/kT)  .  On  the 
basis  of  this  one  would  expect  the  ratio  of  the  mean  height  at  20°C  to  that  at  130°C  would  be 
about  (403/293)^^1 .9.  The  experimentally  determined  ratios  for  the  liquids  tha*  we  have  con¬ 
sidered  fall  into  the  range  of  about  1.0  to  1.5. 
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The  Mean  Height  to  Which  the  Droplets  Rise  Multiplied  hv  Density  to  the  5/6 
Power  Plotted  Versus  the  Square  of  the  Bulk  Modulus.  The  Solid  Line  Is  Drawn 
According  to  Equation  D7).  The  data  shown  are  for  the  six  liquids  listed  in  table  1 
( where  the  various  svtnbols  are  identified)  at  a  temperature  of  20°C. 


Figure  3. 


MODULUS/GPa 


0  12  3 


-  2 

[modulus  /GPaj 


The  Mean  Height  to  Which  the  Droplets  Rise  Multiplied  hv  Density  to  the  5/6 
Power  Plotted  Versus  the  Square  of  the  Bulk  Modulus.  The  Solid  Line  Is  Drawn 
According  to  Equation  D7).  The  data  shown  are  for  the  six  liquids  listed  in  table  1 
( where  the  various  svtnbols  are  identified)  at  a  temperature  of  20°C. 


Figure  3. 


Table.  Data  Used  in  Constructing  Figure  3.  All  Data  Are  Given  for  a  Temperature  of  20°C. 
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This  cannot  be  accounted  for  by  considering  the  change  with  temperature  of  the  density  and  the 
modulus.  It  probably  reflects  the  increasingly  important  role  of  the  viscoelastic  responses,  for 
which  the  relaxation  times  will  be  significantly  shorter  at  the  higher  temperature. 


6.  THE  EFFECT  OF  POLYMERIC  ADDITIVES 

Gerber  and  Stuempfle*  also  report  measurements  at  20°C  on  two  samples  of 
tetrachloroethvlene  to  one  of  which  was  added  1  %  polvisobutylene  (PIB),  while  the  other 
contained  1%  polystyrene  (PS).  In  the  former  case  <z>was  increased  by  nearly  a  factor  of  four 
while  in  the  latter  case  the  increase  was  less  than  5 0% .  Phrased  another  way,  the  effect  of  PIB 
is  nearly  to  double  the  average  droplet  radius,  whereas  PS  produces  only  about  a  20%  increase  in 
droplet  size.  The  model  of  comminution  that  we  have  considered  would  lead  one  to  expect  an 
increase  in  droplet  size  on  the  order  of  10%  or  less.  Clearly  some  other  mechanism  is  operative 
in  determining  the  particle  size.  We  suggest  that  the  effect  of  the  polymeric  molecules  is  to 
act  as  agents  that  terminate  the  flaws  hv  inhibiting  their  growth  under  the  applied  stresses. 
The  PIB  molecules,  which  were  of  substantially  greater  average  chain  length,  are  much  more 
effective  in  this  regard  than  the  (relatively  shorter)  PS  molecules,  and  thus  lead  to  the 
formation  of  rather  large  droplets. 


7.  SUMMARY  AND  CONCLUSIONS 

In  this  report  we  have  presented  a  physical  model  of  the  explosive  comminution  process 
in  liquids.  The  key  feature  of  the  model  involves  the  development  of  flaws  in  the  liquid  under 
the  action  of  large  explosively  generated  stresses.  It  is  at  these  flaws  that  the  material  fails 
via  a  process  that  we  have  termed  "liquid  fracture."  Guided  hv  previous  molecular  dynamics 
(MD)  studies  that  we  have  carried  out  on  the  nonlinear  response  character  of  liquids  (a  partial 
summary  of  these  is  given  in  the  Appendix),  we  have  proposed  a  highly  simplified  view*  of  this 
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process.  Even  this  simplified  picture  is  relatively  successful  in  yielding  semiquantitative  agree¬ 
ment  with  experimental  data. 

In  order  to  develop  a  truly  quantitative  picture,  it  is  our  belief  that  the  dynamical 
aspects  of  the  liquid  failure  mechanisms  must  be  more  fully  understood.  A  thoroughgoing 
molecular  dynamics  investigation  of  the  time-dependent  structural  response  of  liquids  to  high- 
amplitude  shear  and  tensile  stresses  would  prove  most  useful  in  elucidating  the  essential 
features  of  these  nonlinear  viscoelastic  dynamics  that  operate  in  explosive  comminution.  We 
are  at  present  examining  the  possibilities  of  initiating  such  an  effort. 
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APPENDIX 


MOLECULAR  DYNAMICS 


Not  only  in  the  process  of  comminution  but  also  in  a  rather  wide  variety  of  technological 
applications,  for  example  elastohydrodynamic  lubrication,  shock  loading,  solid  propellant 
failures,  etc.,  liquids  and  amorphous  solids  are  subjected  to  extremely  large  shearing  forces. 
Because  the  microscopic  mechanisms  of  importance  are  not  well  understood,  solutions  to 
problems  involving  nonlinear  behavior  are  quite  limited. 

For  example  in  order  to  predict  shear  failure  (or  tearing)  of  a  liquid  under  the  forces  of 
an  explosion  or  shock  wave,  one  must  be  able  to  compute  the  time  evolution  of  the  stress  in  a 
fluid  element  when  it  encounters  a  simultaneously  applied  time  varying  shear  rate  and  a  time 
varying  pressure.  The  response  of  the  system  under  these  conditions  can,  for  low  shearing  rates, 
be  parameterized  in  terms  of  the  shear  viscosity,  0,  a  shear  rigidity  modulus,  G^,,  and  a 
distribution  of  relaxation  times.  At  the  rather  large  shearing  rates  generally  encountered,  this 
relatively  simple  description  must  be  modified  to  incorporate  the  nonlinear  character  of  the 
dynamical  stress  response.  Some  authors*  suppose  that  this  can  be  done  by  appending  the  idea 
of  a  limiting  shear  stress  to  the  low  amplitude,  i.e.,  linear,  shear  response.  In  order  to  test  this 
hypothesis  and,  more  generally,  to  provide  a  framework  on  which  to  base  analytical  repre¬ 
sentations  of  nonlinear  viscoelastic  behavior,  we  have  used  Molecular  Dynamics  (MD)  to  inves¬ 
tigate  a  model  liquid  system  under  conditions  approximating  those  found  in  at  high  shear  rates. 
Even  though  the  model  system  is  particularly  simple,  we  have  been  able  to  reproduce  the  major 
features  of  shear  response  of  real  liquids.  Because  an  MD  experiment  provides  one  with  a 
complete  microscopic  record  of  the  system’s  evolution  in  time,  one  can  perform  "experiments" 
and  "measure"  aspects  of  its  behavior  that  are  not  accessible  in  conventional  experiments,  but 
which  provide  useful  insights  into  the  origin  of  the  properties  of  liquids  under  extreme  stress 
conditions. 

The  MD  calculation  procedure  observes  a  portion  of  liquid  as  it  is  sheared.  Although  a 
simultaneously  applied  time  dependent  pressure  and  shear  rate  would  more  closely  represent  the 
conditions  of  a  liquid  undergoing  explosive  comminution,  we  have  confined  our  attention  up  to 
now  to  considering  the  response  of  a  model  fluid  to  a  time  varying  shear  rate  at  a  series  of 
densities  (pressures). 

THE  MODEL 

MD  is  a  computer  simulation  technique  for  probing  the  microscopic  behavior  of  a  system 
that  allows  one  to  account  for  the  many-body  nature  of  the  molecular  interactions.  It  involves 
following  some  representative  number  of  molecules,  N,  by  solving  the  classical  equations  of 
motion  on  a  high  speed  computer  to  obtain  their  trajectories  in  time  steps  of  typically  0.01  ps. 
Since  only  rather  a  limited  number  of  molecules  is  treated  (generally  less  than  a  thousand),  the 
molecules  are  confined  in  a  cubic  box  which  is  surrounded  by  images  of  itself  to  avoid  severe 
boundary  effects.  In  the  studies  described  below  the  artificialities  introduced  by  this  procedure 
are  rather  insignificant.  From  the  results  of  an  MD  "experiment"  —  essentially  a  record  of  the 
system's  path  through  phase  space  —  one  can  determine  any  measurable  property  of  the  system, 
as  well  as  many  that  are  experimentally  inaccessible. 


*Bair,  S..  and  Winer,  W.  O.  a  Rheological  Model  for  Elastohydrodynamic  Contacts  Based  on 
Primary  Laboratory  Data,  ASME  J.  Lub.  Tech.,  101,  258  (1979). 
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The  model  system  we  consider  here  is  an  assembly  of  108  Lennard-Jones  (LJ) 
particles;  these  molecules  interact  through  pairwise  additive  forces  described  by  the  LJ  6-12 
potential, 


^(r)  =  -4u  (a/r)°  -  (a/r’ 


i2J. 


(1) 


where  r  is  the  separation  of  an  interacting  pair,  and  u  and  a  are  constants  with  dimensions  of 
energy  and  distance,  respectively.  All  the  computed  quantities  are  in  so-called  LJ  reduced  units 
which  are  given  in  terms  of  u,  a,  and  m,  the  mass  of  an  individual  particle.  Although  no  real 
liquid  can  be  characterized  by  a  potential  of  the  LJ  form  (only  the  noble  gas  fluids  Ar,  Kr,  etc., 
are  quantitatively  represented  in  this  way)  their  structural  and  dynamical  behavior  can  be 
expected  to  be  (and  we  have  found  actually  are)  qualitatively  similiar  to  that  of  the  model 
fluid.  As  a  result  one  can  expect  that  the  microscopic  insights  derived  from  MD  experiments 
will  serve  as  useful  guidelines  for  understanding  and  modeling  the  behavior  of  real  liquids  and 
perhaps  for  the  "molecular  engineering"  of  improved  comminution  properties  of  liquids.  For  the 
purpose  of  comparison  with  experiment  frequent  conversion  from  reduced  to  "real"  units  will  be 
made,  using  the  parameters  for  Ar  which  are  given  in  table  A-l. 


Table  A-l.  A  Summary  of  the  Reduced  Units  in  Terms  of  Which  All  Quantities  Are  Given. 

Boltzmann’s  constant  is  denoted  by  kg. 


Ouantity 

Reduced  Unit 

SI  Unit  for  Ar 

distance 

a 

0.3405  nm 

mass 

m 

6.64  x  kg 

energy 

u 

1.65  x  10"21  J 

time 

a(m/u)3^2 

2.16  ps 

densitv 

a’3 

42.1  kg  mol/m3 

temperature 

U/kB 

119.8  K 

pressure,  stress 

u/a3 

41.8  MPa 

modulus  ! 

u/a3 

41.8  MPa 

viscosity 

(mu^/a2 

9.03  x  10~5  Pa  s 

thermal  conductivity 

kp  (m/u)  J/^a"2 

1.88  x  10“2  J  m-1  K-1  s_1 

Two  methods  of  shearing  the  box  of  molecules  cam  be  used.  The  results  of  each  metho/1 
agree  within  statistical  error;  however,  each  method  is  well  suited  to  examine  a  particular 
aspect  of  the  investigated  phenomenon. 


The  first  method  attempts  to  mimic  a  situation  in  which  the  fluid  forms  a  thin  film. 
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agree  within  statistical  error;  however,  each  method  is  well  suited  to  examine  a  particular 
aspect  of  the  investigated  phenomenon. 

The  first  method  attempts  to  mimic  a  situation  in  which  the  fluid  forms  a  thin  film. 
This  thin  film  is  achieved  in  the  model  system  by  employing  periodicity  in  the  x  and  y  directions 
onlv;  the  external  shear  forces  are  applied  by  tvro  fluid  layers  translating  with  the  desired  wall 
velocities,  Uj  and  U^.  Molecules  from  the  three  regions  (i.e.,  the  two  wall  regions  and  the 
liquid  itself)  are  not  allowed  to  mix  and  are  kept  from  doing  so  by  reflection  boundary 
conditions  in  the  z  direction.  The  fluid  walls  (FW)  are  maintained  at  a  preset  temperature  bv 
adjusting  their  velocities  at  each  time  step  of  the  integration  scheme. 
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The  second  method  is  a  modification  of  the  Homogeneous  Shear  technique*  (MHS),  which 
is  better  suited  to  study  the  viscoelastic  behavior  of  the  fluid  because  it  employs  periodicity  in 
all  directions  so  that  a  bulk  material  is  sheared.  It  also  enables  a  desired  strain  rate  to  be 
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MD  ’  M  K A S U K  E  V  E  NT"  OF  SHEAR  VISCOSITY  OF  HIGH  STRAIN  RATES 
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large  values  of  e  were  chosen  since  it  is  desirable  to  explore  the  system’s  behavior  at  shear 
rates  on  the  order  of  the  reciprocal  of  the  viscoelastic  relaxation  times.  For  the  LJ  fluids  these 
are  approximately  unity  (about  1C  ^  s).  The  shear  stress  needed  to  define  a  shear  viscosity  is 
readily  derived  from  the  velocities  and  positions  of  the  molecules.  "Die  a8  component  of  the 
stress  tensor  is  defined  below7, 
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where  V  is  the  volume  of  the  MD  cell,  v?  is  the  o  velocity  component  of  molecule  i  relative  to 
the  average  velocitv  and  or  -  is  the  a  component  of  the  separation  between  molecules  i  and  j. 
The  subscripts  for  the  shear  stress  component,  ^X2>  are  dropped  in  further  discussion. 

At  low  shear  rates  plots  of  shear  stress  versus  strain  rate,  *  ,  are  linear;  however,  at 
high  stress  levels  ( ~  G  / 1 00 )  the  shear  stress  increases  less  rapidly  than  «  .  In  other  words  the 
shear  viscositv  (=•  a/*  )  decreases  from  the  equilibrium  (or  zero  shear  rate)  value  of  r?0  with 
increasing  shear  rate.  The  MHS  method  was  used  to  obtain  tj  versus  i  for  a  variety  of 
densities.  The  forms  of  V/Vq  versus  stress  in  figure  A-2  resemble  closelv  those  obtained  from 
twin-disk  experiments  for  the  fluids  C>P4E,  and  the  lubricants  L63/1271  and  Oxilube  85/140  after 
making  the  appropriate  conversions  to  real  units.  The  ratio  V /V0  descends  markedly  from  unity 
in  the  region  —  50MPa.  The  decrease  in  rj  with  c  is  not  due  to  thermal  heating  because  the 
calculations  w^re  conduoted  isothermallv. 


*  He  yes,  D.  M.  et  a!.,  op.  cit. 
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I  igure  A-l.  The  2  Dependence  of  the  Average  x  Velocity  for  the  p  =  0.8442  and  t  =  0.337  3 
State  Using  the  FW  Model.  The  z  coordinate  system  is  measured  from  the 
center  of  the  film. 
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Figure  A-Z. 


The  Nonlinear  Behavior  of:  A  BP  L63/1271  (P  =  0.8  GPa,  T  =  27UC,  U  =  1.12 
m/s),  O  5P4E  (P  =  0.45  GPa,  T  =  27°C,  U  =  0.6  m/s)  and  rOxiluhe  85/140  (P  = 
1.2  GPa,  T  =30°C,  U  =  2.2  m/s)  from  Twin-Disk  Experiments.  The  MHS 
Model:  ▲  p  =  0.8442,  ▼  p  -  0.92862,  •  p  =  1.01304  and  open  square  p  = 
1.2663. 
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CALCULATION  OF  SHEAR  INDUCED  STRUCTURAL  CHANGES 

The  results  in  the  previous  section  suggest  that  the  model  and  real  liquids  might  share  a 
<  mechanism  of  "failure"  at  high  levels  of  shear  stress.  With  the  support  of  the  previous 

results  we  now  consider  associated  structural  changes  produced  by  the  large  shear  rates  in  the 
lluids  in  order  to  better  understand  the  failure  mechanism, 

A  directional  number  density  function,  n(a),  has  been  calculated  for  the  purpose  of 
investigating  structural  anisotropy.  The  function  n(x),  for  example,  is  the  average  number  of 
mob -rales  from  a  reference  particle  in  the  distance  range  |x|  to  |x|  +  Ax,  where  Ax  was  taken  as 
(  '200,  selected  from  the  N-l  particles  in  the  re-orientated  MD  cell.  Figure  A-3  shows  the  n(a) 
f  the  unsheared  liquid  subtracted  from  that  of  the  sheared  system,  5n(a).  Although  the  maxi- 
•M'un  displacement  is  10%  at  most,  there  is  evidence  of  a  tendency  of  the  molecules  to  re-order 
! lit  msel ves  intv)  xv  layers  (layers  normal  to  the  velocity  gradient).  The  xz  shearing  increases 
tm*  probability  of  the  molecules  in  an  xz  plane  being  found  at  multiples  of  the  intermolecular 
ammeter  from  t-avl,  other  in  the  z  direction.  It  is  energetically  favorable  for  these  layers  to  he 
staggered  in  the  y  direction.  Although  not  so  well  defined,  these  trends  are  present  at  the 
lower  densities.  The  (grosslv  simplified)  picture  that  emerges  is  one  of  corrugated  sheets  of 
atonic  sliding  past  one  another.  A  pictorial  representation  of  this  suggested  structure  is  given  in 
figure  A  1.  Th  »  structural  reorganization  results  from  the  action  of  the  imposed  velocity 
profile  wiiich,  at  the  shear  rates  considered,  dominates  the  molecular  motion.  Major  move¬ 
ments  art-  confined  to  the  line  of  shear. 


Another  aspect  of  the  shear  induced  structural  distortion  is  presented  by  the  time 
dependent  pair  radial  distribution  function,  g  (r,t),  which  is  conveniently  decomposed  into  the 
following  component  configurational  averages: 

R  O  >  t )  =  g(j(r)  +  (  x2  /  r  2  )  uxx  (  r  ,  t )  +  ( v2 / r 2 ) uyy ( r  ,  t  ) 

+  ( z2 / r2 ) uzz ( r , t )  +  ( xz / r2 ) uxz ( r  ,  t ) 


=  . . . ,  (3) 

where  gQ<r)  is  the  pair  radial  distribution  function  of  the  unsheared  medium. 

2 

if  (r,t)  is  the  average  of  (  °ij  ^  the  rac^a*  element  r_^r  +  dr  about  a 

molecule  i  at  time  t,  then 

=  15  foo2r,t*  '  *■  p  r2f1r)’ 


and 


ga/3(r,t)  =  15  f  (r,t)  /  (47rpr2dr),  (4) 

•  here  a  =£8 

The  fa$(r,t)  have  been  calculated  at  various  times  during  the  segment.  The  difference 
between  the  ga,*(r,t)  of  the  steady  state  sheared  system  and  the  unsheared  svstem  at  the  same 
density,  6garr'r,t),  is  shown  in  figure  A-5.  It  reveals  that  shearing  produces  a  net  movement  of 
particles  in  each  coordination  shell  towards  the  origin  molecule. 

The  -.hear  rate  dependent  gxz(r>t)  at  steady  state  of  the  p  =  0.84  svstem  are  compared 
with  the  spherically  averaged  g(r,t)  from  these  samples  in  figure  A-6.  Thev  give  convincing 
vidence  of  angular  distortion  in  each  coordination  shell  such  that  in  the  positive  xz  quadrants 


Appendix 


22 


CALCULATION  OF  SHEAR  INDUCED  STRUCTURAL  CHANGES 

The  results  in  the  previous  section  suggest  that  the  model  and  real  liquids  might  share  a 
<  mechanism  of  "failure"  at  high  levels  of  shear  stress.  With  the  support  of  the  previous 

results  we  now  consider  associated  structural  changes  produced  by  the  large  shear  rates  in  the 
lluids  in  order  to  better  understand  the  failure  mechanism, 

A  directional  number  density  function,  n(a),  has  been  calculated  for  the  purpose  of 
investigating  structural  anisotropy.  The  function  n(x),  for  example,  is  the  average  number  of 
mob -rales  from  a  reference  particle  in  the  distance  range  |x|  to  |x|  +  Ax,  where  Ax  was  taken  as 
(  '200,  selected  from  the  N-l  particles  in  the  re-orientated  MD  cell.  Figure  A-3  shows  the  n(a) 
f  the  unsheared  liquid  subtracted  from  that  of  the  sheared  system,  5n(a).  Although  the  maxi- 
•M'un  displacement  is  10%  at  most,  there  is  evidence  of  a  tendency  of  the  molecules  to  re-order 
! lit  msel ves  intv)  xv  layers  (layers  normal  to  the  velocity  gradient).  The  xz  shearing  increases 
tm*  probability  of  the  molecules  in  an  xz  plane  being  found  at  multiples  of  the  intermolecular 
ammeter  from  t-avl,  other  in  the  z  direction.  It  is  energetically  favorable  for  these  layers  to  he 
staggered  in  the  y  direction.  Although  not  so  well  defined,  these  trends  are  present  at  the 
lower  densities.  The  (grosslv  simplified)  picture  that  emerges  is  one  of  corrugated  sheets  of 
atonic  sliding  past  one  another.  A  pictorial  representation  of  this  suggested  structure  is  given  in 
figure  A  1.  Th  »  structural  reorganization  results  from  the  action  of  the  imposed  velocity 
profile  wiiich,  at  the  shear  rates  considered,  dominates  the  molecular  motion.  Major  move¬ 
ments  art-  confined  to  the  line  of  shear. 


Another  aspect  of  the  shear  induced  structural  distortion  is  presented  by  the  time 
dependent  pair  radial  distribution  function,  g  (r,t),  which  is  conveniently  decomposed  into  the 
following  component  configurational  averages: 

R  O  >  t )  =  g(j(r)  +  (  x2  /  r  2  )  uxx  (  r  ,  t )  +  ( v2 / r 2 ) uyy ( r  ,  t  ) 

+  ( z2 / r2 ) uzz ( r , t )  +  ( xz / r2 ) uxz ( r  ,  t ) 


=  . . . ,  (3) 

where  gQ<r)  is  the  pair  radial  distribution  function  of  the  unsheared  medium. 

2 

if  (r,t)  is  the  average  of  (  °ij  ^  the  rac^a*  element  r_^r  +  dr  about  a 

molecule  i  at  time  t,  then 

=  15  foo2r,t*  '  *■  p  r2f1r)’ 


and 


ga/3(r,t)  =  15  f  (r,t)  /  (47rpr2dr),  (4) 

•  here  a  =£8 

The  fa$(r,t)  have  been  calculated  at  various  times  during  the  segment.  The  difference 
between  the  ga,*(r,t)  of  the  steady  state  sheared  system  and  the  unsheared  svstem  at  the  same 
density,  6garr'r,t),  is  shown  in  figure  A-5.  It  reveals  that  shearing  produces  a  net  movement  of 
particles  in  each  coordination  shell  towards  the  origin  molecule. 

The  -.hear  rate  dependent  gxz(r>t)  at  steady  state  of  the  p  =  0.84  svstem  are  compared 
with  the  spherically  averaged  g(r,t)  from  these  samples  in  figure  A-6.  Thev  give  convincing 
vidence  of  angular  distortion  in  each  coordination  shell  such  that  in  the  positive  xz  quadrants 


Appendix 


22 


Figure  A-4.  \  Pictorial  Representation  of  the  Structural  Changes  That  Take  Place  on  Going 

from  an  Unsheared  (a)  to  a  Sheared  (b)  State,  The  arrows  denote  the  line  and 
plane  of  shear. 
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the  i dM tie  of  en<  h  shell  is  on  average  depleted  of  particles  when  compared  with  the  unshoared 
liquid.  In  contrast,  the  outside  of  each  shell  has  an  excess  of  molecules  in  the  positive  xz  quad- 
when  compared  with  the  unsheared  medium.  The  opposite  changes  take  place  in  the 

ti ve  quadrants. 


FFI  i’(  i  OF  STRUCTURAL  REORGANIZATION  ON  MOLECULAR  MOBILITY 


U  is  to  he  expected  that  the  effect  of  the  shear  should  not  be  confined  to  structural 
’  h  liges  hut  must  also  alter  the  self-dif  f  usional  dynamics  of  the  molecules.  To  assess  this,  the 
x,  v,  and  z  components  of  the  mean  square  displacements,  <(r  (t))>,  which  exclude  the  imposed 

flow,  have  been  calculated,  i.e.,  T 
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lection  of  trio  states  are  presented  table  A-2.  The  value  of  Da  ~  0.032  for  the  unsheared 
li  pm.  near  the  triple  point  compares  favorably  with  the  equilibrium  value  of  0.033.  The  Da  i.f 
the  unsheared  liquid  decrease  by  approximately  an  order  of  magnitude  on  increasing  the  density 
tro.n  p  0.8442  to  1.2nt>3.  However,  at  larger  shear  rates  the  density  has  a  diminishing  in  flu- 
'*nce  on  the  diffusion  coefficients.  The  high  shear  stresses  enhance  the  fluiditv  and  create  a 
o  ju id  which  has  directionality.  All  the  diffusion  coefficient  components  rise  with  shear  rate 
at  a  particular  density  because  the  structural  reorganization  that  takes  place  on  shearing  appar¬ 
ently  c  rentes  paths  along  which  particles  can  more  readily  move.  The  D  are  tvpicallv  two- 
thirds  <  f  the  Dx  and  02,  indicating  that  the  self-diffusion  is  favored  in  the  shearing  plane. 
1  bus,  paradoxic  ally,  although  structural  aspects  of  the  the  liquids  under  shear  present  a  more 
solid-like  appearance,  other  more  dynamically  related  properties  such  as  self-diffusion  and 
shear  viscosity  manifest  changes  which  are  associated  with  enhanced  fluiditv. 

Thus  we  see  that  the  high  stresses  cause  extreme  distortions  of  the  system's  structure. 
I  ne  inter  molecular  forces  are  not  sufficiently  strong  to  support  the  enormous  stress  levels  that 
develop;  consequently  there  is  a  rupturing  and  reorganization  of  the  structure  resulting  in 
enhanced  flow  behavior,  due  to  the  alignment  of  molecules  into  "glide’  planes  along  the  line1  of 
shear. 


These  structural  rearrangements  are  also  manifest  in  changes  in  the  so-called  normal 
pressure  components,  PUCl  (-  -  Oaot)  Tvpicallv,  Pzz  >  Pxx  >  P^r  although  their  differences  are 
small  when  compared  with  each  component's  change  on  shearing.  In  figure  A- 7  the  time  depen- 
dem  e  of  the  increase  in  (t)  is  shown.  It  is  found  that  the  time  scale  of  the  development  ol 

the  normal  pressure  components  is  longer  than  that  of  the1  shear  stress. 

Hie  inference  in  the  viscoelastic  and  structural  reorganization  times  can  he  better 
umh  cstood  b\  further  consideration  of  the  radial  distribution  functions. 

1  he  tine  dept  ndence  of  goto,  and  g^p  (r,t)  reveal  that  structural  reorganization 

saving  the  symmetry  of  shear,  ^z/r  >,  is  much  faster  than  that  of  <otd/r^>  symmetry.  Within 
;h  sern  ,,  and  <z^/r^>  evolve  more  rapidlv  to  their  steady  state  values 

fh  o.  <yu/r^  For  example,  when  P  =  1.01304  and  Gxz  -  0.7108,  the  former  three  averages  have 
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umh  cstood  b\  further  consideration  of  the  radial  distribution  functions. 

1  he  tine  dept  ndence  of  goto,  and  g^p  (r,t)  reveal  that  structural  reorganization 

saving  the  symmetry  of  shear,  ^z/r  >,  is  much  faster  than  that  of  <otd/r^>  symmetry.  Within 
;h  sern  ,,  and  <z^/r^>  evolve  more  rapidlv  to  their  steady  state  values 

fh  o.  <yu/r^  For  example,  when  P  =  1.01304  and  Gxz  -  0.7108,  the  former  three  averages  have 
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nearly  reached  their  steady  state  values  by  t  =  0.25  whereas  that  of  the  <y  /r  >  summation 
hardly  differs  from  zero  by  that  time.  Thus  we  find  that  when  a  liquid  is  sheared,  the  shearing 
forces  alter  its  structure  so  that  the  steady  state  stress  attained  corresponds  to  a  new  lower 
viscosity,  not  that  characterizing  the  starting  liquid.  At  high  rates  of  shear  the  structural 
evolution  necessary  to  attain  this  new  viscosity  is  significantly  slower  than  the  viscoelastic 
shear  relaxation.  Viewed  simply  this  means  that  the  original  viscosity  can  remain  for  a  suffi¬ 
ciently  long  time  to  enable  the  shear  stress  to  exceed  its  steady  state  value. 


Table  A  2.  The  Density  and  Shear  Rate  Dependence  of  the  Directional  Diffusion  Coefficients, 

Da  ,  Obtained  from  the  MD  Experiments. 
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it  has  come  to  be  increasingly  recognized  that  the  process  of  fracture  involves  time 
dependent  nonlinear  viscoelastic  behavior.  In  order  to  better  understand  the  time  dependent 
nonlinear  viscoelasticity,  the  MD  calculations  were  conducted  with  a  time  varying  strain  rate. 
In  what  follows  we  describe  the  viscoelastic  MD  “experiments"  and  the  results.  What  is  most 
encouraging  is  that  in  every  case  the  behavior  of  real  liquids  is  mimicked  by  the  MD  results. 

We  consider  the  dynamical  response  of  the  LJ  material  to  an  imposed  shear  that  is 
'switched  on"  suddenly  using  the  Homogeneous  Shear  method.  Two  types  of  MD  experiments 
are  1  interest:  (a)  the  response  of  the  system  to  a  steady  shearing  rate  (uniform  velocity 
gradient)  €  -  d  v  /  d  z  that  is  initiated  at  time  t  =  0;  and  (b)  the  response  of  the  system  to  a 
constant  shear  strain  (a  pulse  of  shearing  rate)  imposed  at  time  t  =  t’.  Linear  response  theory 
provides  us  with  relations  between  the  two  types  of  responses:  the  shear  stress  response  to  a 
step  of  shear  strain  can  be  written  as 


<j(t)/e  =  G^cMt),  (6) 

where  (t)  is  the  normalized  (<f>(0)  -  1)  shear  stress  relaxation  function.  In  terms  of  these  pa¬ 
rameters  the  shear  stress  response  to  a  “small”  steady  shear  rate  imposed  at  t  =  0  is, 


a  (t)/e  =  Gc 


dt  ’  <t>  ( t  *  )  . 


(7) 


Since  the  steady  state  (t-  +  oo)  value  of  this  is  just  the  viscosity  of  the  material,  rj  ,  we  have, 

%  -  °-/r dt’  * (r)  =  G<»T- 

The  last  equality  serves  to  define  the  shear  realization  time  r. 
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Experiments  of  the  type  described  as  (a)  above  were  then  carried  out  as  a  function  of 
t  for  values  of  €  well  beyond  the  range  for  which  linear  behavior  can  be  expected.  A  represen¬ 
tative  selection  of  the  results  is  displayed  in  figures  A-8  and  A-9  in  which  the  density  (con¬ 
stant  f_ )  and  shear  rate  (constant  density)  dependence  of  17  are  shown,  respectively. 

Several  features  of  these  figures  are  worthy  of  comment.  At  short  times  the  system 
behaves  largely  as  an  elastic  solid  in  that  the  stress  increases  linearly  with  time, 

0(t)  =  G  ^et.  (9) 


The  shear  rigidity  modulus  rises  linearly  with  pressure  at  roughly  the  same  rate  as  more  compli¬ 
cated  fluids,  as  is  revealed  in  figure  A-10.  Liquids  cannot  permanently  support  a  shear  strain 
and  consequently  the  stress  levels  off  to  a  value  determined  by  the  shear  viscosity  of  the 
sheared  state.  There  is  evident  "shear  thinning,"  or  a  decrease  of  viscosity  with  increasing 
hear  rate.  The  predicted  value  of  the  stress  build-up  for  the  p  =  0.92862,  G  =  0.6963  state 
using  equation  (7),  shown  in  figure  A-8,  clearly  overestimates  the  actual  stress  obtained. 

A  iso,  a  maximum  in  che  stress  versus  time  curve  for  those  simulations  in  the  high 
t».  nsitv/sh  ar  rate  regime  is  observed.  This  indicates  that  the  large  shear  deformations  are 
reducing  structural  distortions  which  lag  in  time  somewhat  behind  the  growth  of  the  shear 
stress  but  follow  rather  closely  the  observed  time  variation  of  the  normal  pressure  components 
in  the  materia!  This  relatively  slower  development  of  what  can  be  thought  of  as  transient 
stress-induced  glide  planes  is  consistent  with  the  observed  shear  thinning  and  stress  "overshoot" 
behavior  presented.  It  is  also  consistent  with  the  observation  that  the  stress  relaxation  function 
decays  more  rapidly  in  a  highly  sheared  system  than  in  an  equilibrium  system  (see  discussion  of 
gt .  ^(t)  above).  A  comparison  of  the  two  types  of  behavior  is  presented  in  figure  A-ll.  The  form 
of  these  i  t)  consists  of  a  rapidly  (t  <  0.2)  decaying  portion  which  is  presumably  inertial  in 
origin  and  due  to  small  motions  of  the  molecules  within  the  random  network  structure.  Further 
stress  relaxation  requires  the  cooperative  rearrangement  of  larger  molecular  groups  and  is 
1  msoquentlv  slower.  The  attenuation  of  4>(t)  with  increasing  shear  rate  is  indicative  of  the 
enhanced  stress  relieving  properties  of  the  structurally  reorganized  system.  The  stress  over¬ 
shoot  observed  in  nonlinear  viscoelastic  experiments  on  polymers  has  a  form  that  is  quite 
similar  to  the  curves  obtained  for  the  LJ  system.  While  the  microscopic  mechanisms  for  the 
two  systems  are  obviously  quite  different,  the  general  principle  that  associates  the  stress  over¬ 
shoot  phenomena  (as  well  as  shear  thinning)  with  a  time  dependent  structural  reorganization 
process  remains  quite  valid. 

In  this  context  it  is  perhaps  relevant  to  note  that  a  more  rapid  decay  of  the  stress 
relaxation  function  can  be  thought  of  as  both  a  suppression  of  the  long  relaxation  times 
associated  with  stress  relaxation  and  also  a  reduction  in  the  shear  rigidity  modulus,  at  low 
frequencies,  G,  (a  rather  well  known  phenomenon  in  polymeric  systems).* 

Thus  we  see  that  the  viscoelastic  studies  at  high  shear  rates  reflect  the  structural  reor¬ 
ganization  that  is  ultimately  responsible  for  failure  or  fracture  of  the  liquid.  The  fact  that 
every  property  "measured"  at  high  shear  rates  are  found  in  real  and  far  more  complicated 
liquids  suggests  strongly  that  MD  studies  can  be  utilized  to  investigate  the  comminution  process 

m  those  liquids. 


*  Y.  (1.  Yanovski  in  Relaxation  Phenomena  in  Polymers,  ed.  G.  Bartenev  and  Y.  V.  Zelenev 
(Khimiva,  Leningrad,  1972),  p.  113. 
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The  Time  Dependence  of  the  Reduced  Stress,  a(t)/e: _ p  =  0.8442,  go  = 

.  (a)  Actual  .  (b)  Predicted  Using  Equation  (4)  p  =  0.92862, 

0.6963; _ p  =  1.01304,  eQ  =  0.7168;  _ P  =  1.2663,  gq  =  0.7721. 
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Figure  A-9. 


Figure  A-9. 


l*he  Pressure  Dependence  of  :  ASantotrac  40,  O  Di(2-ethylhexyl)  Phthalate 
from  a  Twin-Disk  Apparatus.  Solid  line  G^  =  (0.293  +  1.70P)  GPa  di(2- 
e)phthalate.  The  filled  in  O  is  obtained  from  MHS  MD. 
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